# select only one of the least pruned solutions

args = commandArgs(TRUE) 
file1 = args[1] # adjacency matrix of the initial hairball
file2 = args[2] # all best solutions of the contextualization
cutoff = args[3]

G = read.csv(file1,header=TRUE,stringsAsFactors=FALSE,quote="",sep="\t",check.names=FALSE); G = as.matrix(G)
B = read.csv(file2,sep=" ",header=FALSE,stringsAsFactors=FALSE,quote="")

ind_optim = which(G!=0)  # nonzero edges in the raw adjacency matrix (it goes through the matrix in columns)
ind_actinh = which(G==2) # unspecified edges in the raw adjacency matrix
if (length(ind_actinh)+length(ind_optim) != nrow(B)) stop("length(ind_actinh)+length(ind_optim) != nrow(B)")

M = matrix(nrow=length(ind_optim), ncol=1) #[edges,1]
I.G = which(G!=0,arr.ind=TRUE) #coordinates of ind_optim in raw adjacency matrix
rownames(M) = sprintf("%d__%d",I.G[,1],I.G[,2])	# rownumber__colnumber (where rows and cols are genes in the raw adjacency matrix)

I.G = cbind(I.G,'interact'=0)

# select one solution
least_pruned = which.max(apply(B[1:length(ind_optim),,drop=FALSE],2,sum)) #max number of edges retained (ind_optim positions=1)


b_least = t(B[,least_pruned])

extractTint=function(b) {
	G1 = G  
	G1[ind_actinh[!!b[(length(ind_optim)+1):length(b)]]] = 1	#assign signs to unspecified edges
	G1[ind_actinh[!b[(length(ind_optim)+1):length(b)]]] = -1
	G1[ind_optim[!b[1:length(ind_optim)]]] = 0 #remove edges

	I.G[,3] = G1[I.G[,1:2]]

	a1 = sprintf("%d__%d",I.G[,1],I.G[,2])
	rownames(I.G) = a1
	M[,1] = I.G[match(rownames(M), rownames(I.G)),3]

	act = NULL
	inh = NULL
	# 1. activation
	I = which(G1==1,arr.ind=TRUE)
	if(length(I)!=0){
	  act = cbind(colnames(G1)[I[,1]], "-->", colnames(G1)[I[,2]])
	}
	# 2. inhibition
	I = which(G1==-1,arr.ind=TRUE)
	if(length(I)!=0){
	  inh = cbind(colnames(G1)[I[,1]], "--|",colnames(G1)[I[,2]])
	}  

	interactions=rbind(act,inh)
	interactions
}

int_least=extractTint(b_least)

write.table(int_least, file=sprintf("interactions_least_pruned-cutoff%s.txt",cutoff),row.names=FALSE,quote=FALSE, col.names=FALSE, sep="\t")

# check size of the GRN; GRN<10 is not considered enough
grn_size = length(union(int_least[,1],int_least[,3]))
print(sprintf("The contextualized GRN contains %s TFs",grn_size))
if (grn_size < 10) print("NOT BIG ENOUGH")
